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ABSTRACT 

The extended gamma ray source MGRO J 1908+06, discovered by the Milagro air shower detector in 2007, 
has been observed for ^4 years by the ARGO-YBJ experiment at TeV energies, with a statistical signifi- 
cance of 6.2 standard deviations. The peak of the signal is found at a position consistent with the pulsar PSR 
J1907+0602. Parametrizing the source shape with a two-dimensional Gauss function we estimate an extension 
o'ext = 0.49°±0.22°, consistent with a previous measurement by the Cherenkov Array H.E.S.S.. The observed 
energy spectrum is dN/dE = 6.1±1.4 x 10~ 13 (E/4 TeV)" 2 ' 54±a36 photons cm" 2 s" 1 TeV -1 , in the energy range 
^1-20 TeV. The measured gamma ray flux is consistent with the results of the Milagro detector, but is ^2-3 
times larger than the flux previously derived by H.E.S.S. at energies of a few TeV. The continuity of the Mi- 
lagro and ARGO-YBJ observations and the stable excess rate observed by ARGO-YBJ along 4 years of data 
taking support the identification of MGRO J 1908+06 as the steady powerful TeV pulsar wind nebula of PSR 
J 1907+0602, with an integrated luminosity above 1 TeV ~1.8 times the Crab Nebula luminosity. 
Subject headings: gamma rays: general - pulsars: individual (MGRO J1908+06) 



1. INTRODUCTION 

The Galactic gamma ray source MGRO J1908+06 was dis- 
covered by the Milagro air shower detector in a survey of the 
Galactic plane at a median energy of ~ 20 TeV (Q]). The data 
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were consistent both with a point source and with an extended 
source of diameter <2.6°. Assuming a spectrum ocE" 2 - 3 , the 
measured flux at the median energy of 20 TeV is 8.8±2.4 x 
10~ 15 photons cm" 2 s" 1 TeV" 1 . 

A marginal detection of a source consistent with the posi- 
tion of MGRO J 1908+06 was already reported by the Tibet 
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AS-7 array ( 1331) . but not confirmed in a more recent paper 

The source was later observed by the H.E.S.S. ® and 
VERITAS (32) Cherenkov telescopes. In particular, H.E.S.S 
detected an extended source (HESS J1908+063) at energies 
above 300 GeV (~1 1 standard deviations of statistical signif- 
icance) positionally consistent with MGRO J1908+06. The 
measured source extension, evaluated assuming a symmetri- 
cal two-dimensional Gaussian shape, was a ext = 0.34°^q q3. 

H.E.S.S. reported a power law differential energy spectrum 
with a photon index of 2.10 ± 0.07 s(fl , ± 0.2 jys in the energy 
range 0.3-20 TeV, and a flux at 1 TeV of (4.14 ± 0.32 stat ± 
0.83 JW ) x 10" 12 photons cm" 2 s" 1 TeV" 1 . The integrated flux 
above 1 TeV is 17% that of the Crab Nebula. 

After the release of the Bright Source List by the Fermi 
collaboration (3), Milagro reported the association of MGRO 
J1908+06 to the LAT pulsar 0FGL J1907.5+0602 (later re- 
named PSR J1907+0602), pulsating with a period of 106.6 
ms @. The peak of the Milagro emission was 0.3° off the 
pulsar, but consistent with the pulsar location within the mea- 
surement error (0.27°). Assuming a spectrum ocE" 2 6 , Mila- 
gro reported a flux ofll6.7±15.8 x 10" 17 photons cm" 2 s" 1 
TeV" 1 , at the median energy of 35 TeV. 

The association of MGRO J 1908+06 with PSR J 1907+0602 
was also supported in (5), where a multiwavelength study of 
the pulsar and the surrounding region has been performed 
with radio, X-ray and Fermi gamma ray data. Because of the 
small angular distance between the pulsar and the centroid of 
the H.E.S.S. extended source, the authors argue that the latter 
is plausibly the Wind Nebula of the pulsar. 

Performing an off-pulse measurement, Fermi set an upper 
limit to the HESS J1908+063 flux in the energy region 0.1-25 
GeV, suggesting that the spectrum has a low-energy turnover 
between 20 GeV and 300 GeV. With radio and X-ray data, a 
lower limit to the pulsar distance was set to ^3.2 kpc, deriving 
for the nebula a physical size >40 pc. 

Later, Milagro evaluated the energy spectrum of the source 
in the 2-100 TeV region, reportin g a hard p ower law spec- 
trum with an exponential cutoff {Smith 2009). The best fit 
obtained is dN/dE = 0.62 x 10" 11 E" 150 exp(-E/14.1) pho- 
tons cm" 2 s _1 TeV" 1 , where E is the energy in TeV. This flux is 
in disagreement with that given by H.E.S.S. at a level of 2-3 
standard deviations, being about a factor 3 higher at 10 TeV. 
The authors suggest that the discrepancy can be simply due to 
a statistical fluctuation, or to the fact that Milagro, given its 
relatively poor angular resolution, integrates the signal over a 
larger solid angle compared with H.E.S.S., and likely detects 
more of the diffuse lateral tails of the extended source. 

In this work we report on the observation of MGRO 
J 1908+06 with the ARGO-YBJ detector performed during the 
years 2007-201 1. After a brief description of the detector and 
a detailed presentation of the data analysis technique, we re- 
port our results concerning the extension and the energy spec- 
trum of the source. 

2. THE ARGO-YBJ EXPERIMENT 

The ARGO-YBJ detector is located at the Yangbajing Cos- 
mic Ray Laboratory (Tibet, China) at an altitude of 4300 m 
above sea level. It consists of a ~74x 78 m 2 carpet made 
of a single layer of Resistive Plate Chambers (RPCs) with 
^92% of active area, sorrounded by a partially instrumented 
(^20%) area up to ~ 100 x 1 10 m 2 . The apparatus has a modu- 
lar structure, the basic data acquisition element being a cluster 



(5.7x7.6 m 2 ), made of 12 RPCs (2.8x1.25 m 2 ). The RPCs 
are operated in streamer mode by using a gas mixture (Ar 
15%, Isobutane 10%, TetraFluoroEthane 75%) suitable for 
high altitude operation. 

Each RPC is read by 80 strips of 6.75x61.8 cm 2 (the spa- 
tial pixels), logically organized in 10 independent pads of 
55 .6 x 6 1 .8 cm 2 which are individually acquired and represent 
the time pixels of the detector dToh . In addition, in order to 
extend the dynamical range up to PeV energies, each RPC is 
equipped with two large size pads (139 x 123 cm 2 ) to collect 
the total charge developed by the particle hitting the detector 
d25l) . The full experiment is made of 153 clusters for a total 
active surface of ^6600m 2 . 

ARGO-YBJ operates in two independent acquisition 
modes: the shower mode and the scaler mode (fill) . In this 
analysis we refer to the data recorded from the digital read- 
out in shower mode. In this mode, an electronic logic has 
been implemented to build an inclusive trigger, based on a 
time correlation between the pad signals, depending on their 
relative distances. In this way, all the shower events giving a 
number of fired pads N pa d > N tr j g in the central carpet in a 
time window of 420 ns generate the trigger. This trigger can 
work with high efficiency down to N £I ^=20, keeping negligi- 
ble the rate of random coincidences (Il5l) . 

The time of each fired pad in a window of 2 /isec around the 
trigger time and its location are recorded and used to recon- 
struct the position of the shower core and the arrival direction 
of the primary particle. 

In order to perform the time calibration of the 18360 pads, 
a software procedure has been developed, based on the Char- 
acteristic Plane method (1241) which using the secondary parti- 
cles of large vertical showers as calibration beams, iteratively 
reduces the differences between the measured times and the 
temporal fit of the shower front (1121) . 

The full detector is in stable data taking since 2007 Novem- 
ber with the trigger condition N„-^=20 and a duty cycle ~ 
86%. The trigger rate is ~3.5 kHz with a dead time of 4%. 

3. DETECTOR PERFORMANCE 

The angular resolution and the pointing accuracy of the de- 
tector have been evaluated by using the Moon shadow, i.e. 
the deficit of cosmic rays in the Moon direction, observed by 
ARGO-YBJ with a statistical significance of ^9 standard de- 
viations per month. The shape of the shadow provides a mea- 
surement of the detector Point Spread Function (PSF), and its 
position allows the individuation of possible pointing biases. 

The data have been compared with the results of a Monte 
Carlo simulation which describes the propagation of cosmic 
rays in the Earth magnetic fields, the shower development in 
the atmosphere by using the CORSIKA code (1231) . and the 
detector response with a code based on the GEANT package 
J22h . The PSF measured with the cosmic rays has been found 
in excellent agreement with the Monte Carlo evaluation, con- 
firming the reliability of the simulation procedure ( 18). 

The angular resolution for gamma rays is evaluated by sim- 
ulating the events from a gamma ray source with a given spec- 
trum and daily path in the sky. It results smaller by ^30-40% 
compared with the angular resolution for cosmic rays, due 
to the better defined time profile of the showers. In general, 
the PSF for gamma rays can be described by the sum of two 
Gaussian distributions. For a Crab-like source, the radius of 
the opening angle which optimizes the signal-to-background 
ratio for events with N pad >60 (300) is 0.86° (0.44°) and con- 
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FIG. 1 . — PSF-smoothed significance map of the MGRO J1908+06 region 
obtained by ARGO-YBJ, for events with N pai j >60. Open circle: position of 
the center of MGRO J 1908+06, as measured by Milagro. The error bars give 
the linear sum of the statistical and systematic errors. Open triangle: centroid 
of HESS J1908+063. Black filled circle: Fermi pulsar PSR J1907+0602. 
White filled circles: Fermi gamma ray sources, according to the 2nd Fermi 
Catalogue 1 27). White crosses: TeV sources detected by H.E.S.S. in the same 
region. The white line represents the Galactic plane. 



tains ~50% of the signal. 

The Moon Shadow has also been used to check the abso- 
lute energy calibration of the detector, by studying the west- 
ward shift of the shadow due to the geomagnetic field. The 
observed displacement as a function of the event multiplicity 
Npad is in excellent agreement with the results of the Monte 
Carlo simulation. From this analysis the total absolute en- 
ergy scale error, including systematic effects, is estimated to 
be less than 13% (El). 

4. DATA ANALYSIS AND RESULTS 

At the ARGO-YBJ site MGRO J 1908+06 culminates at the 
zenith angle of 24° and is visible for 5.38 hours per day with 
a zenith angle less than 45°. The dataset used in this analysis 
refers to the period from November 2007 to December 201 1 
and contains all the showers with zenith angle less than 45° 
and Npad >20. The total on-source time is 6867 hours. 

To study the gamma ray emission from a source, a 16° x 16° 
sky map in celestial coordinates (right ascension and declina- 
tion) with 0. 1 ° x 0. 1 ° bin size, centered on the source position, 
is filled with the detected events. 

In order to extract the excess of gamma rays, the cosmic ray 
background has to be estimated and subtracted. 

The time swapping method (1141) is used to evaluate the 
background: for each detected event, n "fake" events (with 
n = 10) are generated by replacing the original arrival time 
with new ones, randomly selected from an event buffer which 
spans a time T of data taking. Changing the time, the fake 
events maintain the same declination of the original event, but 
have a different right ascension. With these events a new sky 
map (background map) is built, with a statistics n times larger 
than the "true" event map in order to reduce the fluctuations. 
To avoid the inclusion of the source events in the background 



evaluation, the showers inside a circular region around the 
source (with a radius related to the PSF and depending on 
Nparf) are excluded from the time swapping procedure. A cor- 
rection of the number of swaps is made to take into account 
the rejected events in the source region (1211) . The value of the 
swapping time T is ~ 3 hours, in order to minimize the sys- 
tematic effects due to the environmental parameter variations. 

In order to extract the source signal the maps are smoothed 
according to the detector PSF, determined by Monte Carlo 
simulations for different N pa d intervals. Finally, the smoothed 
background map is subtracted to the smoothed event map, ob- 
taining the "excess map", where for every bin the statistical 
significance S of the excess is given by: 

S=(N „-N of f)/y/6Nt n + 5N2 ff 

with N on = £j Ni Wi and N ff = W; In. In these ex- 

pressions Ni and B,- are the number of events of the i' h bin of 
the "event map" and "background map", respectively, w, is a 
weight proportional to the value of the PSF at the angular dis- 
tance of the \ th bin, and n is the number of swaps. The sum 
is over all the bins inside a radius R, chosen to contain the 
PSF. Since the number of events per bin is large, the fluctu- 
ations follow the Gaussian statistics, hence the errors on N on 
and N ff are: 5N „ = ^/T^w] and SN /f = ^/£;B;W? /n 2 . 

In order to study the signal in different energy regions, 
the maps are built for 8 different N pac { intervals, namely 20- 
39, 40-59, 60-99, 100-199, 200-299, 300-499, 500-999, and 
>1000. These maps are then combined to have "integral 
maps" for different N pal j thresholds. 

Analysing the data recorded in 4 years, the sky maps of 
the MGRO J 1908+06 region show a significant excess at the 
source position for different N pa d thresholds. The larger sig- 
nificance is given by events with N pal j > 20, with 7.3 standard 
deviations. Increasing N pa d, the significance decreases. For 
Npa,/ > 1000 no signal is present. 

The distributions of the significances outside the source re- 
gion follow a standard normal distribution, showing the cor- 
rectness of the background evaluation procedure. 

As will be discussed in Section 4.3, the signal with 
N pa rf=20-59 is largely affected by the Galactic diffuse gamma 
ray flux, and only events with N pa d > 60 will be used in the 
study of the source morphology and flux. Fig.l shows the 
significance map for events with N pa d > 60, where the source 
signal reaches 6.2 standard deviations. 

Studying the source on the time scale of one year, the an- 
nual excess rate results to be consistent with the total av- 
erage rate, indicating that the gamma ray flux from MGRO 
J 1908+06 is likely due to a steady emission. 

4. 1 . Source position and extension 

To evaluate the position and extension of the source, the 
events with N pac / > 60 are used. We assume a source shape 
described by a symmetrical two-dimensional Gaussian func- 
tion with r.m.s. a exr . Fitting the non-smoothed excess map 
to a function given by the convolution of the above Gaussian 
and the detector PSF, we found the best-fit position at RA. = 
19&08T and decl. = 6°24\ with a statistical error of 12' and 
a systematic error of 6' per axis. The position found is consis- 
tent with the the Milagro measurement and with the centroid 
of HESS J1908+063 (RA. = 19 A 07 m 54 s and decl. = 6°16'7", 
with a statistical error of 2.4' and a systematic error of 20" per 
axis). 



4 



Bartoli et al. 



1 14000 
> 

^12000 
o 




1 1.5 2 2.5 

Angular distance (deg) 

FIG. 2. — Number of excess events with N pflI ; >60 as a function of the an- 
gular distance from the best-fit centroid position, compared with the expected 
distributions for different source extensions. 



The value of a ext that best fits the data is 0.49°±0.22°, con- 
sistent with the H.E.S.S. estimation of 0.34°. 

Fig. 2 shows the distribution of the angular distance from 
the best-fit centroid position, compared with the simulated 
distributions corresponding to the extensions a ext =0A9° and 
ov.if=0°. The two curves are normalized to the same number 
of excess events. 

4.2. Energy spectrum 

In order to study the energy dependence of the signal the 
events are divided in different subsets according to the number 
of hit pads, and a sky map is built for each subset. For this 
analysis we define 4 intervals: N pad = 20-59, 60-199, 200- 
499 and N pa d >500. The intervals have been chosen in order 
to have a signal of comparable statistical significance. 

For the spectrum evaluation we assume a power law depen- 
dence: dN/dE = K E" 7 . The values of K and 7 are derived by 
comparing the number of the excess events detected in each 
of the previously defined N pac { intervals with the correspond- 
ing ones given by simulations assuming a set of test spectra. 
The reliability of this procedure has been tested studying the 
the Crab Nebula signal (fl3h . 

For each N pa d interval, the number of excess events is ob- 
tained integrating the sky map around the source position up 
to a distance ip max , where ipmax is the radius of the opening 
angle which maximizes the signal to background ratio. The 
value of ip max depends on the source extension a ex , and on the 
detector PSF, and is provided by simulations. For the exten- 
sion we use the value a exl = 0.49°, according to our measure- 
ment. On the other hand, the PSF for a given N /M ^ interval 
is not precisely determined, since it depends both on the de- 
tector characteristics and on the spectrum index 7, which is 
unknown. To solve this "circular" problem an iterative proce- 
dure has been applied. 

First, an initial index 7 = 2.5 is assumed, and the corre- 
sponding values of ip max for every N pa d interval are deter- 
mined via simulations. The number of events observed in 
ipmax are then used to evaluate a new spectral slope 7 which 
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FIG. 3. — Gamma ray flux from MGRO J 1908+06 measured by different 
detectors. ARGO-YBJ: the continuous line (colored red in the online version) 
is the best fit to data. The dashed area represents the one standard deviation 
error. H.E.S.S.: the dotted line is the best fit to the points 1:9). Milagro 1: 
flux value assuming a spectrum ocE -2,3 Q]). Milagro 2: flux value assuming 
a spectrum <xE~ 2,6 @). Milagro 3: the dashed line (colored blue in the online 
version) is the spectrum fit according to Smith (2009) and the vertical lines 
are the errors (at one standard deviation) for some values of the energy. The 
plotted errors are purely statistical for all the detectors. 



is returned to the first step, to calculate a new set of ipmax, and 
so on. Given the relatively weak dependence of the PSF on 
7, a small number of iterations is sufficient to terminate the 
process successfully and provide the parameters of the best fit 
spectrum. 

4.3. Contribution from diffuse flux 

Since the source is located on the Galactic plane, the ob- 
served flux could be affected by the diffuse gamma ray emis- 
sion produced by cosmic rays interacting with the matter and 
the radiation fields of the Galaxy. Given the relatively large 
opening angles used in the measurement, the photons from 
the diffuse radiation falling in the observational window of 
MGRO J 1908+06 could artificially increase the flux detected 
from the source direction. The amount of this contribution can 
be evaluated by analysing the data collected from the Galactic 
plane region close to the source. 

The flux of very high energy Galactic gamma rays in the 
region of MGRO J1908+06 (Galactic coordinates I = 40.39° 
and b = -0.79°) is poorly known. 

The first evidence of a diffuse Galactic emission at TeV 
energies has been reported by the Milagro detector (1171) . A 
significant dependence of the flux on the Galactic latitude 
and longitude has been found in a later analysis with events 
of median energy 15 TeV (2). The same paper reports the 
expected energy spectrum for two different sectors of the 
Galactic plane, for energies from 10 keV to 100 TeV, ac- 
cording to the GALPROP model (HH HI), "optimized" to fit 
the measurements by EGRET in the 40 MeV - 10 GeV en- 
ergy range and by Milagro at 15 TeV. Concerning the MGRO 
J 1908+06 region, the expected average flux at 1 TeV in the 
area of Galactic coordinates / e [30°,65°] and b € [-2° ,2°] is 
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~2x 10~ 9 photons TeV -1 cm" 2 s" 1 sr" 1 . 

A preliminary flux measurement at energies E >300 GeV, 
obtained with the ARGO-YBJ data, is reported by Ma 
(201 1), who derives the average gamma ray spectrum for / £ 
[25°,65°] and b £ [-2°,2°]. This estimate results lower but 
still consistent with the expectation of the above model. 

For our purposes, given the variation of the emission along 
the Galactic plane and its strong dependence on the latitude, it 
is preferable to evaluate the diffuse flux in a restricted region, 
adjacent to the source position. 

We consider two sky regions, L; and L2, of size Al = 5° 
and Ab = 2 x ip max (where ip max depends on the N pac i interval) 
whose centers have the same latitude of MGRO J 1 908+06 and 
are located at both sides of the source, at a longitudinal dis- 
tance of 5.5°. 

Analysing the ARGO-YBJ data, a global excess of statis- 
tical significance 3.0 standard deviations over the cosmic ray 
background is observed in Li + L2 for events with N pat j >20. 
This excess is interpreted as due to the diffuse Galactic emis- 
sion plus the contribution of 5 gamma ray sources discovered 
by H.E.S.S., namely HESS J1912+101 (i), W49B ©, HESS 
J1857+026 and HESS J1858+020 0), and HESS J1849-000 

The individual fluxes of these objects are below the ARGO- 
YBJ sensitivity, while the total flux is (29±3)% that of the 
Crab Nebula at 1 TeV (ED). In particular, HESS J1912+101 
and HESS J1857+026 have a flux -10% and -17% that of 
the Crab Nebula, respectively. The number of events from 
these sources expected to fall in Li and L2 is evaluated via 
simulations, using the fluxes measured by H.E.S.S., and gives 
a global contribution of (40±14)% to the observed diffuse ex- 
cess. After the subtraction of this contribution, taking into 
account the different exposures of Li and L2, we evaluate the 
number of events due to the diffuse emission expected to fall 
into the observational window of MGRO J1908+06. 

We found that the ratio between the number of events ex- 
pected from the diffuse emission and those observed from 
the source direction is Rj = 0.33±0.18 for showers with 
N pa rf=20-59, and R4 <0.15 (at one sigma level) for showers 
with Npad >60. 

As a comparison, the values derived by using the "opti- 
mized" GALPROP diffuse emission model given by Abdo et 
al.(2008) are R d = 0.57 for N pfld =20-59, and R d = 0.23 for 
Nporf >60. These values are larger than those obtained with 
the ARGO-YBJ data. It should be noted, however, that the 
above model is based on a measurement by Milagro which 
does not take into account all the gamma ray sources located 
in the studied region, and could overestimate the diffuse flux. 

The larger contribution of the diffuse emission for N pa d=2Q- 
59 is due to the wider opening angle used in this interval 
(VWt=2.0 o ). Because of these estimates, to avoid a possible 
large systematic effect in the flux evaluation, we restrict our 
spectral analysis to the events with N pa d >60. 

Performing the procedure described in Section 4.2, we 
fit the data of the 3 intervals N /Mrf =60-199, 200-499 and 
Nporf >500. The best fit spectrum obtained is: dN/dE = 
6.1±1.4 x 10" 13 (E/4 TeV)- 2 ' 54±0 ' 36 photons cm" 2 s" 1 TeV -1 , 
valid in the energy region 1-20 TeV. The median energies cor- 
responding to the 3 N pa d intervals are 2.4, 5.1 and 12.8 TeV, 
respectively. 

As a comparison, if we do not exclude the data with 
^^=20-59, the best-fit spectrum is: dN/dE = 1.36±0.29 
x 10" 12 (E/3 TeV)- 2 ' 65±0 - 25 photons cm" 2 s" 1 TeV" 1 , which 



gives a flux 2 1 % higher at E = 1 TeV. 

Beside the statistical errors and the systematics due to the 
diffuse contribution discussed before, our measurement could 
be affected by an additional systematic error mainly due to 
the background evaluation, to the absolute energy scale deter- 
mination, to the pointing accuracy, to environmental effects 
and to the Monte Carlo simulations, for a global effect that 
we estimate to be < 30% (Aielli et al. 2010). 

In the case of an extended source, a possible further cause 
of systematics could be the uncertainty in the extension, and 
the consequent use of an incorrect opening angle to extract 
the signal. Therefore we have also evaluated the spectrum 
assuming a ex , = 0.34°, as measured by H.E.S.S. The resulting 
flux differs from the previous one by less than 5% in the whole 
energy range considered in the analysis. 

The obtained spectrum is shown in Fig. 3, together with 
those reported by H.E.S.S. and Milagro. The flux is signif- 
icantly higher than that given by H.E.S.S. in the 1-10 TeV 
energy range , but is consistent with the Milagro spectrum 
(Smith 2009). The hard spectrum with exponential cutoff ob- 
tained by Milagro produces a worse, but still acceptable, fit to 
our data. Given the reduced significance of the excess at high 
energies, we are not able to constrain the shape of the spec- 
trum above 10 TeV and to definitively rule out a high energy 
cutoff. 

5. DISCUSSION AND CONCLUSIONS 

The gamma ray source MGRO J 1908+06 has been stud- 
ied by ARGO-YBJ analyzing —4 years of data. An excess 
with significance 6.2 standard deviations is observed in a po- 
sition consistent with previous measurements by Milagro and 
H.E.S.S.. The peak of the signal occurs at R.A. = \9 h 0% m \ s 
and decl. = 6°24' (with statistical and systematic errors of 
~0.2° and 0.1° per axis, respectively) and lies at a distance of 
22' from PSR J 1907+0602, consistent with the pulsar location 
within the measurement error. 

The signal is due to emission from an extended region. Af- 
ter taking into account the detector PSF, the extension of the 
source is found to be a ext = 0.49°±0.22°. 

The photon spectrum in the range 1-20 TeV follows a sim- 
ple power law with a spectral index 2.54±0.36, though a 
harder spectrum with a high energy cutoff cannot be ruled 
out. 

T he spectrum i s found to be consistent with the Milagro re- 
sult (ISmifh 2009 ) but not with the H.E.S.S. best fit in the 1-10 
energy range, the flux measured by ARGO-YBJ at 4 TeV be- 
ing a factor 2.6 larger. At - 20 TeV the ARGO-YBJ, H.E.S.S. 
and Milagro fluxes are consistent within the errors, and are 
also in agreement with the first Milagro measurement (Qjl. 

Since a contribution to this measurement is expected from 
the Galactic diffuse emission, data from two sky regions lo- 
cated at both sides of the source and centered at the same lati- 
tude have been used to determine this contamination. Accord- 
ing to this estimate, the diffuse Galactic gamma ray emission 
is expected to contribute to the signal above 1 TeV for less 
than — 15%, and cannot account for the observed disagree- 
ment. 

Being the difference with H.E.S.S. at the level of 2.5 stan- 
dard deviations, the discrepancy could be simply due to statis- 
tical fluctuations, or to the combination of statistical and sys- 
tematic uncertainties. However these latter have been accu- 
rately studied, giving a global error on the flux less than 30%. 
Indeed, the spectrum of the Crab Nebula obtained by ARGO- 
YBJ results to be in good agreement with the Cherenkov de- 
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tectors measurements (Q [3lb . The extension of the source 
should not give such an additional systematic error to explain 
the observed difference. 

On the other hand, a similar discrepancy is found in the 
observation of the extended source MGRO J2031+41, located 
in the Cygnus region, for which ARGO-YBJ (fl9h and Milagro 
© report a flux significantly larger than that measured by the 
Cherenkov Telescopes MAGIC and HEGRA. 

In principle, one cannot exclude the possibility of a flux 
variation as the origin of the observed disagreement among 
the detectors. Milagro, H.E.S.S. and ARGO-YBJ data have 
been recorded in different periods. Milagro integrates over 
seven years (July 2000 - November 2007) while the total 
H.E.S.S. data set only amounts to 27 hours of sparse obser- 
vations during 2005-2007, before the ARGO-YBJ measure- 
ment. However, a possible flux variation seems unlikely, since 
the average fluxes measured by Milagro and ARGO-YBJ in 
two contiguous periods covering a total time of 1 1 years, are 
consistent. 

Moreover, it should be noted that if MGRO J 1908+06 is 
the pulsar wind nebula associated to PSR J 1907+0602, the 
gamma ray emission originates from a region whose size has 
been estimated to be >40 pc (5), implying that the variation 



time scale cannot be less than ~130 years, unless relativistic 
beaming effects are present. 

In conclusion, MGRO J 1908+06 is observed by ARGO- 
YBJ as a stable extended source, likely the TeV nebula of 
PSR J 1907+0602, with a flux at 1 TeV -67% that of the Crab 
Nebula. Assuming a distance of 3.2 kpc, the integrated lu- 
minosity above 1 TeV is —1.8 times that of the Crab Nebula, 
making MGRO J 1908+06 one of the most luminous Galactic 
gamma ray sources at TeV energies. 
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